Prepare Tidy CVC Datasets (with SWMM model hydrology)

Setup the basic working environment


In [ ]:
%matplotlib inline

import os
import sys
import datetime
import warnings
import csv

import numpy as np
import matplotlib.pyplot as plt
import pandas
import seaborn
seaborn.set(style='ticks', context='paper')

import wqio
import pybmpdb
import pynsqd

import pycvc

min_precip = 1.9999
big_storm_date = datetime.date(2013, 7, 8)

palette = seaborn.color_palette('deep', n_colors=6)
pybmpdb.setMPLStyle()
POCs = [
    p['cvcname'] 
    for p in filter(
        lambda p: p['include'],
        pycvc.info.POC_dicts
    )
]

if wqio.testing.checkdep_tex() is None:
    tex_msg = ("LaTeX not found on system path. You will "
               "not be able to compile ISRs to PDF files")
    warnings.warn(tex_msg, UserWarning)
    
warning_filter = "ignore" 
warnings.simplefilter(warning_filter)

Load External Data (this takes a while)


In [ ]:
bmpdb = pycvc.external.bmpdb(palette[3], 'D')
nsqdata = pycvc.external.nsqd(palette[2], 'd')

Load CVC Database


In [ ]:
cvcdbfile = "C:/users/phobson/Desktop/scratch/cvc/cvc.accdb"
cvcdb = pycvc.Database(cvcdbfile, nsqdata, bmpdb)

Hydrologic Relationships

$V_{\mathrm{runoff, \ LV1}} = \max\left(0,\: -12.05 + 2.873\, D_{\mathrm{precip}} + 0.863 \, \Delta t \right)$


In [ ]:
def LV1_runoff(row):
    return max(0, -12.0 + 2.87 * row['total_precip_depth'] + 0.863 * row['duration_hours'])

ED-1

$\log \left(V_{\mathrm{runoff, \ ED1}}\right) = 1.58 + 0.000667 \, I_{\mathrm{max}} + 0.0169 \, D_{\mathrm{precip}} $

$V_{\mathrm{bypass, \ ED1}} = \max \left(0,\: -26.4 + 0.184 \, I_{\mathrm{max}} + 1.22 \, D_{\mathrm{precip}} \right)$

$V_{\mathrm{inflow, \ ED1}} = \max \left(0,\: V_{\mathrm{runoff, \ ED1}} - V_{\mathrm{bypass, \ ED1}} \right)$


In [ ]:
def ED1_runoff(row):
    return 10**(1.58 + 0.000667 * row['peak_precip_intensity'] + 0.0169 * row['total_precip_depth'] )

def ED1_bypass(row):
    return max(0, -26.4 + 0.184 * row['peak_precip_intensity'] + 1.22 * row['total_precip_depth'])

def ED1_inflow(row):
    return max(0, ED1_runoff(row) - ED1_bypass(row))

LV-2

$\log \left(V_{\mathrm{runoff, \ LV2}}\right) = 1.217 + 0.00622 \, I_{\mathrm{max}} + 0.0244 \, D_{\mathrm{precip}} $

$V_{\mathrm{bypass, \ LV2}} = 0$

$V_{\mathrm{inflow, \ LV2}} = \max \left(0,\: V_{\mathrm{runoff, \ LV2}} - V_{\mathrm{bypass, \ LV2}} \right)$


In [ ]:
def LV2_runoff(row):
    return 10**(1.22 + 0.00622 * row['peak_precip_intensity'] + 0.0244 * row['total_precip_depth'] )

def LV2_bypass(row):
    return 0

def LV2_inflow(row):
    return max(0, LV2_runoff(row) - LV2_bypass(row))

LV-4

$\log \left(V_{\mathrm{runoff, \ LV4}}\right) = 1.35 + 0.00650 \, I_{\mathrm{max}} + 0.00940 \, D_{\mathrm{precip}} $

$V_{\mathrm{bypass, \ LV4}} = \max \left(0,\: 7.37 + 0.0370 \, I_{\mathrm{max}} + 0.112 \, D_{\mathrm{precip}} \right)$

$V_{\mathrm{inflow, \ LV4}} = \max \left(0,\: V_{\mathrm{runoff, \ LV4}} - V_{\mathrm{bypass, \ LV4}} \right)$


In [ ]:
def LV4_runoff(row):
    return 10**(1.35 + 0.00650 * row['peak_precip_intensity'] + 0.00940 * row['total_precip_depth'] )

def LV4_bypass(row):
    return max(0, 7.36 + 0.0370 * row['peak_precip_intensity'] + 0.112 * row['total_precip_depth'])

def LV4_inflow(row):
    return max(0, LV4_runoff(row) - LV4_bypass(row))

Water quality loading relationship

$ M_{\mathrm{runoff}} = V_{\mathrm{runoff}} \times \hat{\mathbb{C}}_{\mathrm{inflow}}\left(\mathrm{landuse,\ season}\right) $

$ M_{\mathrm{bypass}} = V_{\mathrm{bypass}} \times \hat{\mathbb{C}}_{\mathrm{inflow}}\left(\mathrm{landuse,\ season}\right) $

$ M_{\mathrm{inflow}} = M_{\mathrm{runoff}} - M_{\mathrm{bypass}} $

$ M_{\mathrm{outflow}} = V_{\mathrm{outflow}} \times \mathbb{C}_{\mathrm{outflow}} $

Define the site object for the reference site and compute its median values ("influent" to other sites)


In [ ]:
LV1 = pycvc.Site(db=cvcdb, siteid='LV-1', raingauge='LV-1', tocentry='Lakeview Control', 
                 isreference=True,  minprecip=min_precip, color=palette[1], marker='s')

LV1.runoff_fxn = LV1_runoff

Lakeview BMP sites get their "influent" data from LV-1


In [ ]:
def rename_influent_cols(col):
    if col.lower() in ['parameter', 'units', 'season']:
        newcol = col.lower()
    else:
        newcol = 'influent {}'.format(col.lower())
        
    return newcol.replace(' nsqd ', ' ').replace(' effluent ', ' ')

LV_Influent = (
    LV1.medians("concentration", groupby_col='season')
       .rename(columns={'effluent stat': 'median'})
       .rename(columns=rename_influent_cols)
)

LV1.influentmedians = LV_Influent
LV_Influent.head()

Elm Drive's "influent" data come from NSQD


In [ ]:
ED_Influent = (
    cvcdb.nsqdata
         .seasonal_medians
         .rename(columns=rename_influent_cols)
)
ED_Influent.head()

Remaining site objects


In [ ]:
ED1 = pycvc.Site(db=cvcdb, siteid='ED-1', raingauge='ED-1',
                 tocentry='Elm Drive', influentmedians=ED_Influent, 
                 minprecip=min_precip, isreference=False,
                 color=palette[0], marker='o')
ED1.runoff_fxn = ED1_runoff
ED1.inflow_fxn = ED1_inflow


LV2 = pycvc.Site(db=cvcdb, siteid='LV-2', raingauge='LV-1',
                 tocentry='Lakeview Grass Swale', influentmedians=LV_Influent, 
                 minprecip=min_precip, isreference=False,
                 color=palette[4], marker='^')
LV2.runoff_fxn = LV2_runoff
LV2.inflow_fxn = LV2_inflow

LV4 = pycvc.Site(db=cvcdb, siteid='LV-4', raingauge='LV-1',
                 tocentry=r'Lakeview Bioswale 1$^{\mathrm{st}}$ South Side', 
                 influentmedians=LV_Influent, 
                 minprecip=min_precip, isreference=False,
                 color=palette[5], marker='v')
LV4.runoff_fxn = LV4_runoff
LV4.inflow_fxn = LV4_inflow

Fix ED-1 storm that had two composite samples


In [ ]:
ED1.hydrodata.data.loc['2012-08-10 23:50:00':'2012-08-11 05:20', 'storm'] = 0
ED1.hydrodata.data.loc['2012-08-11 05:30':, 'storm'] += 1

Replace total inflow volume with estimate from simple method for 2013-07-08 storm


In [ ]:
storm_date = datetime.date(2013, 7, 8)
for site in [ED1, LV1, LV2, LV4]:
    bigstorm = site.storm_info.loc[site.storm_info.start_date.dt.date == storm_date].index[0]
    inflow = site.drainagearea.simple_method(site.storm_info.loc[bigstorm, 'total_precip_depth'])
    site.storm_info.loc[bigstorm, 'inflow_m3'] = inflow
    site.storm_info.loc[bigstorm, 'runoff_m3'] = np.nan
    site.storm_info.loc[bigstorm, 'bypass_m3'] = np.nan

Export project-wide tidy datasets

Hydrologic (storm) data

The big event from July 8, 2013 is retained in this step


In [ ]:
hydro = pycvc.summary.collect_tidy_data(
    [ED1, LV1, LV2, LV4], 
    lambda s: s.tidy_hydro
).pipe(pycvc.summary.classify_storms, 'total_precip_depth')
hydro.to_csv('output/tidy/hydro_swmm.csv', index=False)

Water quality data

The loads from the big event on July 8, 2013 are removed in this step


In [ ]:
wq = (
    pycvc.summary
        .collect_tidy_data([ED1, LV1, LV2, LV4], lambda s: s.tidy_wq)
        .pipe(pycvc.summary.classify_storms, 'total_precip_depth')
        .pipe(pycvc.summary.remove_load_data_from_storms, [big_storm_date], 'start_date')
)
wq.to_csv('output/tidy/wq_swmm.csv', index=False)

Individual Storm Reports

(requires $\LaTeX$)


In [ ]:
for site in [ED1, LV1, LV2, LV4]:
    print('\n----Compiling ISR for {0}----'.format(site.siteid))
    site.allISRs('composite', version='draft')